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Abstract 

We investigate numerically a models of the static spherically symmetric boson-fermion stars 
in scalar-tensor theory of gravity with massive dilaton field. The proper mathematical model of 
such stars is interpreted as a nonlinear two-parametric eigenvalue problem with unknown internal 
boundary. We employ the Continuous Analogue of Newton Method (CANM) which leads on 
each iteration to two separate linear boundary value problems with different dimensions inside 
and outside the star, respectively. Along with them a nonlinear algebraic system for the spectral 
parameters - radius of the star Rs and quantity Q is solved also. 

In this way we obtain the behaviour of the basic geometric quantities and functions describing 
dilaton field and matter fields which build the star. 
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1 Introduction 

In this paper we present a numerical method for solving the equations of the general scalar-tensor 
theories of gravity including a dilaton potential term for the general case of mixed boson-fermion star. 
This method is an impruvment of the method concerning the similar problem proposed in our recent 
work The original domain is splitted to two domains: inner inside the star and external outside the 
star. The solutions in these regions are obtained separately and after that they are matched. Before 
we go to the substantial part of this paper we will describe briefly the physical problem. 

2 Main model 

Boson stars are gravitationally bound macroscopic quantum states made up of scalar bosons jl^, ^ 
1^, ^]. They differ from the usual fermionic stars in that they are only prevented from collapsing 
gravitationally by the Heisenberg uncertainty principle. For self-interacting boson field the mass of the 
boson star, even for small values of the coupling constant turns out to be of the order of Chandrasekhar's 
mass when the boson mass is similar to a proton mass. Thus, the boson stars arise as possible 
candidates for non-baryonic dark matter in the universe and consequently as a possible solution of 
the one of outstanding problems in today's astrophysics - the problem of nonluminous matter in the 
universe. Most of the stars are of primordial origin being formed from an original gas of fermions and 
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bosons in the early universe. That is why it should be expected that most stars are a mixture of both, 
fermions and bosons in different proportions. 

Boson-fermion stars are also good model to learn more about the nature of strong gravitational 
fields not only in general relativity but also in the other theories of gravity. 

The most natural and promising generalizations of general relativity are the scalar-tensor theories 
of gravity ||] . In these theories the gravity is mediated not only by a tensor field (the metric of 

space-time) but also by a scalar field (the dilaton) . These dilatonic theories of gravity contain arbitrary 
functions of the scalar field that determine the gravitational "constant" as a dynamical variable and 
the strength of the coupling between the scalar field and matter. It should be stressed that specific 
scalar-tensor theories of gravity arise naturally as low energy limit of string theory |]l5[ ^ |2^, |oj 
which is the most promising modern model of unification of all fundamental physical interactions. 

Boson stars in scalar-tensor theories of gravity with massless dilaton have been widely investigated 
recently both numerically and analytically j25[ |l6|, |2^, ^ ^, ||, ^ . Mixed boson-fermion stars in scalar 
tensor theories of gravity however have not been investigated so far in contrast to general relativistic 
case where boson-fermion stars have been investigated jl^ . 

In present paper we consider boson-fermion stars in the most general scalar-tensor theory of gravity 
with massive dilaton. 

In the Einstein frame the field equations in presence of fermion and boson matter are: 

V„.V'> + -U'{v) = -^a{ip) (t + T 

~ (1) 



where Vm is the Levi-Civita connection with respect to the metric gmn, ("i = 0, ...,3; n — Q, ...,3). 
The constant is given by — S-kG^ where is the bare Newtonian gravitational constant. The 
physical gravitational "constant" is G*^^((p) where A{(f) is a function of the dilaton field ip depending 
on the concrete scalar-tensor theory of gravity. T4^(^'+^') is the potential of boson field. The dilaton 
potential U{ip) can be written in the form U{lp) = m|jy((p) where m/^ is the dilaton mass and V{if) 
is dimcnsionless function of (p. 

The function a{ip) = ^ [ln^((p)] determines the strength of the coupling between the dilaton field 

B F 

Lp and the matter. The functions T and T are correspondingly the trace of the energy-momentum 

m ^ • ^ . . . 

tensor of the fermionic matter[J and bosonic matter . Their explicit forms are 
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-A\^) \di^+d'^ - 2A\^)W{^+^)\ C , (2) 

(£+p)w,„ti"-pC- (3) 

Here is a complex scalar field describing the bosonic matter while 4"+ is its complex conjugated 
function. The energy density and the pressure of the fermionic fluid in the Einstein frame are e = 
A'^{ip)e and p — A'^{ip)p where e and p are the physical energy density and pressure. Instead to give 
the equation of state of the fermionic matter in the form p = p{e) it is more convenient to write it in 
a parametric form 

s^SogitA P^sofip) (4) 
where Eq is a properly chosen dimensional constant and /i is dimcnsionless Fermi momentum. 

^In the present article we consider fermionic matter only in macroscopic approximation, i.e., after averaging quantum 
fluctuations of the corresponding fermion fields. Thus we actually consider standard classical relativistic matter. 
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The physical four- velocity of the fluid is denoted by u^. The potential for the boson field has the 
form 



2 4 ^ ^ 



Field equations together with the Bianchi identities lead to the local conservation law of the energy- 
momentum of matter 



(5) 



We will consider a static and spherically symmetric mixed boson-fermion star in asymptotic flat 
space-time. This means that the metric gmn has the form 



rfs^ = exp [i/{n)] dt^ - exp [\{n)] dU^ - TZ^ {dO'^ + sin^ 6 dip^) 



(6) 



where TZ, 9, ip are usual spherical coordinates. The field configuration is static when the boson field ^ 
satisfies the relationship 

* = a{TZ) exp(ia;f). 

Here w is a real number and d'{TZ) is a real function. Taking into account the assumption that have been 
made the system of the field equation is reduced to a system of ordinary differential equations. Before 
we explicitly write the system we are going to introduce a rcscalcd (dimensionless) radial coordinate 
by r = rriBTi- where ms is the mass of the bosons. From now on, a prime will denote a differentiation 
with respect to the dimensionless radial coordinate r. After introducing the dimensionless quantities 



CO I— 

0= , u = \/K^a, A = 

rriB 



A 
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and defining the dimensionless energy- momentum tensors as = ^^T^ the components of the 
dimensionless energy-momentum tensor of the fermionic and bosonic matter become correspondingly 



rpO 



B 



F 



Tl=-bA\^)f{n), 



da 



nV^(r)exp(-t.(r)) + — exp(-A(r)) 



dr 



Q^a-^{r)exp{-i/{r)) + 



da_ 
dr 



exp(-A(r)) 



- A\^)W{a% 
- A^{^p)W{(j''). 



Here the parameter b = describes the relation between the Compton length of dilaton and the 
usual radius of neutron star in general relativity. 
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The functions T and T describing the trace of energy-momentum tensor have a form: 



T= -A\v) 



f2V^(r)exp(-!/(r)) - 



fda 



\ dr 



exp(-A(r)) 



- AA^{ip)W{a'^ 



T= bA\^) [.9(m)-3/(m)]. 

For the independent (dimensionless) radial coordinate r we have r G [0, -R^] U [Rs, oo) where < 
i?s < oo is the unknown radius of the fermionic part of the mixed boson-fermion star. 

With all definitions we have given the main system of differential equations governing the structure 
of a static and spherically symmetric boson-fermion star can be written in the following form: 

1. In the interior of the fermionic part of the star (the functions in this domain are subscribed by i) 
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1 - exp(A) 1 ^ 
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--exp(A) 

2 da 
r dr 



dr 



dr 



-crexp(A) exp(-i/) + 2A2(<^)T4^ (cr^) 
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g(M) + /(m) 



-F2 + a{if!) 



dip 
dr 



dip 
dr 



dip 
Ik 



(7) 



Here \(r)^v{r),Lp{r),a{r) and /i(r) are unknown functions of r and is a unknown parameter. 
Having in mind the physical assumptions we have to solve the equations (|^) under following boundary 
conditions: 



A(0) = J(0) = ^(0) = 0, f7(0)=a„ MO)=Mc 



(8) 



KRs) - (9) 

where Uc and Hc are the values of density of the bosonic and fermionic matter, respectively at the star's 
centre. The first three conditions in (^) guarantee the nonsingularity of the metrics and the functions 
A(r), <p(r), a{r) at the star's centre. 

2. In the external domain (subscribed by e) there is not fermionic matter, i.e., one can suppose 
formally that the function /i(r) = if a; > -R^. The fermionic part of the energy-momentum tensor 
vanishes identically also and thus the differential equations with respect to the rest four unknown 
functions A(r), J^(r), ^p{r) and a{r) are: 
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dr 
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l-exp(A) J 
F2,e = r < cxp(A) 



dip 
dr 



r dr 2 dr 



■ exp(A) 



2 da 
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(10) 



-crexp(A) n^exp{-iy) + 2A'^{ip)w'{a^) . 
As it is required by the asymptotic flatness of space-time the boundary conditions at the infinity 



are 



u{co) — 0, ip{ao) — 0, a{co) — 
where we denote (•)(cx3) = limr^oo(' )('')■ 



(11) 
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We seek for a solution [A(r), i^{r), (p{r), a{r),ji{r), Rg , ft] subjected to the nonlinear ODEs (|^) and 
([l0|), satisfying the boundary conditions (^), @ and ([ll]). At that we assume the function /i(r) is 
continuous in the interval [0, whilst the functions A(r), v^r) are continuous and the functions 
(p{r), cr(r) are smooth in the whole interval [0, cx)), including the unknown inner boundary r = i?^. 

The so posed BVP is a two-parametric eigenvalue problem with respect to the quantities Rg and 

n. 

Let us emphasize that a number of methods for solving the free-boundary problems are considered 
in detail in |2|. 

Here we aim at applying the new solving method to the above formulated problem. This method 
differs from that one proposed in |^ and for the governing field equations written in the forms (|^) and 
([To|) it possesses certain advantage. 



3 Method of solution 

At first we scale the variable r using the Landau transformation and in this way we obtain a fixed 
computational domain. Namely 

a; = ^,a; e [0, 1] U [l,oo). 
Rs 

For our further considerations it is convenient to present the systems (|^) and (^o|) in following equivalent 
form as systems of first order ODEs: 

-y[ + RsF,{RgX, y,,n) = 0, (12) 
-y'^ + RsFe{RsX,ye,n) = (13) 

with respect to the unknown vector functions 

yi{x) = {X{x),i^{x),ip{x),£,{x),(T{x),'i]{x),fi{x))'^ , 

ye{x) = {X{x),iy{x),(p{x),£,{x),a{x),T]{x))'^ 

and right hand sides = (Fi, i^2, C, ^3, ^4, ^s)"^, = {Fi, F2,t F^,!], Fi)'^ where (.)' stands for 
differentiation towards the new variable x. 

For given values of the parameters Rs and the independent solving of the inner system ( p!2| ) 
requires seven boundary conditions. At the same time we have at disposal only six conditions of the 
kind (^) and (^) . In order to complete the problem we set additionally one more parametric condition 
(the value of someone from among the functions A(a;), v{x), ip{x),^{x), a{x) or r]{x)) at the point a; = 1). 
Let us set for example 

'P^{l)^Vs (14) 

where (ps is a parameter. Then the boundary conditions (||), (^) and (Q of the inner BVP can be 
presented in the form: 



Ba^^y^{Q) - Do,, = 0, Bi^,y,(l) - Di,,((p,) = 0. (15) 

Here the matrices Bo.i — diag{l,0, 0, 1, 1, 1, 1), i^o.i ~ diag{0, 0, 0, 0, CTc, 0, fic), Bij = diag{0, 0, 1, 0, 0, 0, 1), 
Z?i^, = dm.g(0, 0,(^^0,0,0, 0). 

Obviously the solution in the inner domain x G [0, 1] depends not only on the variable x, but it is 
a function of the three parameters Rs, ips as well, i.e., — yi{x, Q, Rs, tps)- 

In the external domain x > 1 the vector of solutions 

ye{x) = {X{x),i'{x),(p{x),£,{x),a{x),r]{x))'^ 

is 6D. Thereupon six boundary conditions are indispensable to solving of the equation (p^). At the 
same time the three boundary conditions (O) are only known. Let us consider that the solution (x) 
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in the inner domain x G [0, 1] is knowledged. Then we postulate the rest three deficient conditions 
to be the continuity conditions at the point x = l.The first of them is similar to the condition (^) 
and the else two we assign to arbitrary two functions from among A(a;), 1^(2;), £,{x), a{x) and ri{x), for 
example 

Ae(l)-A,(l), (pe(l)=^., CTe(l) =^»(1). 

It is convenient to present the boundary conditions in the external domain in matrix form again: 

Bi^eYeil) - Di^eiVs) ^ 0, Soo,eye(oo) = (16) 

where the matrices Bi^e = diag{l, 0, 1, 0, 1, 0), Di^e = diag{Xi{l), 0, fs, 0, cri(l), 0), -Boo,e = diag{0, 1, 1, 0, 1, 0). 

Let the solutions — yi{x,D,, Rg, tps) and ye — ye{x,fl, Rg, ^Ps) be supposed known. Generally 
speaking for given arbitrary values of the parameters Rg , and ips the continuity conditions with 
respect to the functions i'{x),^{x) and r]{x) at the point x — 1 are not satisfied. We choose the 
parameters Rg,^ and in such manner that the continuity conditions for the functions i'{x), ^(x) 
and rj{x) to be held, i.e., 



Ui,Rg,n,^s)-^,{i,Rs,n,^s) =0, (17) 

These conditions should be interpreted as three nonlinear algebraic equations in regard to the un- 
known quantities Rg , 57 and ips ■ The usual way for solving of the above mentioned kind of equations 
( |l7| ) is by means of various iteration methods, for example Newton's methods. The traditional tech- 
nology similarly to the methods like shutting |Q, requires separate treatment of the BVPs and the 
algebraic continuity equations and brings itself to additional linear ODEs for elements of the corre- 
sponding to ( [itI ) Jacobi matrix. These elements are functions of the variable x and they have to be 
known actually only at the point x = 1. The solving of the both nonlinear BVPs (|l^), ( |l5| ) and (|l3|), 
( |l6|) along with the attached linear equations is another hard enough task. 

At the present work using the CANM |l^ , ||l^ we propose a common treatment of both, differential 
and algebraic problems. 

We suppose that the nonlinear spectral problem (p^, (^Sj), (p^, (16) and ( p7| ) has a "well separated" 
Jilt exact solution. Let the functions yi^o{x),yefl{x) and the parameters Rsfi,i^o,(ps,o are initial 
approximations to this solution. The CANM leads to the following iteration process: 

yi,fc+i(x) =yi.k{x) +TkZi^kix), (18) 

ye,fc+l(2;) =yeM{^) +TkZe,k{x), (19) 

Rs,k+1 = Rs,k+TkPk, (20) 

ftk+1 ^ilk+TkUJk, (21) 

'Ps,k+1 = Vs,k + Tk(j}k- (22) 

Here Tk € (0, 1] is a parameter which can rule the convergence of iteration process. The increments 
2,i^k{x), Ze^k{x), Pkji-^k and 4>k, k = 0,1,2,... satisfy the linear ODEs (for sake of simplicity henceforth 
we will omit the number of iterations k): 

dF / BF \ BF- 

-z: + i?.^z,+ (^i?.^ + F,:jp + i?.^. ^yl~R.F., (23) 

-< + i?.f^z.+ (i?.f|^ + F.)p + i?.^. ^y'^-R.F^. (24) 

All the coefficients and right hand sides as well in the above two equations are known functions of the 
arguments x, Rg, ^ by means of the solution from the previous iteration. We seek for the unknowns 
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2,i{x) of equation (23) and Ze{x) of equation ( p^ as a linear combinations with coefficients p, a; and 0: 

z.i{x) = Si(a::) + pui(a;) + LJVi(a:) + ^Wj(x), (25) 
Ze{x) = Se(a::) + pue(a:;) + LJVe(a;) + 0We(a:). (26) 

Here Si{x),Ui{x), Vi{x), Wi{x), Si{x) and Se(x),Ue{x), Ve(x), We(a;) are new unknown functions, 
which are defined in either, internal or external domains. Substituting for the decomposition ( p5[ ) into 
equation (E3|) after reduction we obtain: 



-u^ + Q,(a;)u, =-fF, + i?,^j 
-w-+Qj(a;)wi =0 



(27) 



where Qi{x) = Rs ^^'(^°^'y"^) stands for a square matrix (7x7), which consists of the Frechet 
derivatives of operator F.^ at the point {yi{x), Rg,^}- 

Similarly applying the CANM to the boundary conditions (^5|) and taking into account the depen- 
dence of matrix Di i on the parameter ips yields: 

BoM^) = Do,^ - So,,y.(0), Bi,,z,(l) = D^^, - Bi,,y,(l) - D[ ,-0. 



By means of the decomposition (|2^) we obtain the following eight boundary conditions (four left + 
four right) for the equations (p7|): 

Bo,zu^(0) = 0, Si,,u,(l)=0, 
Bo,»v,(0) = 0, Si,,v,(l) = 0, 

Bo,^w,{Q) = 0, Si,,w,(l) - -I?'i.,((^s). 



Let us now substitute for decomposition (^6|) into the linear equations for external domain (^) . As 
result we obtain the following four vector equations with regard to the unknown functions Sg (x) , Ug (x) , 
Ve(a;) and We(x) with eight boundary conditions (four left + four right): 

-Se + Qe(a;)se =y'^-RsFe, 



-u[.+Qe{x)Ue ^-{Fe+R. 

-Wg + (5e(a;) We =0 



dF. ''''' (29) 



Here Qe (a;) = Rs ^'^'^^^^^''-^'^'^^ is a square matrix (6 x 6) whose elements are Frechet's derivatives of 
the operator Fe at the point {ye{x),Rs,^}. 

The corresponding linear BC are obtained in the same way as (28) and they become: 

Bi^eSe(l) = i^i^e ~ -Bl,eye(l), Boo,eSe{oo) ^ - Boo^eV e{oo) , 

Bi^eUe(l) =0, Boo,eUe(oo) = 0, 

Bl,eVe(l) =0, Boo,eVe((X3) =0, 

Sl,eWe(l) = -D[,,{lPs), Boo.eWe(oo) = 0. 

In the end to compute the increments {p,uj,4'\ of parameters Rst^ and ips we use the three 
conditions (p7|). 

Let the solutions of linear BVP (p7|), ( psf ) and (|29|), ( ^ ) at the fcth iteration stage are assumed 
to be known. For sake of the simplicity we introduce the vector y(a;) = {v{x),S^{x),rj{x))'^ . For 
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two arbitrary functions hi{x) and he{x), defined in left and right vicinity of the point x = 1, we 
set A/i = /ie(l) — hiC\). Then applying the CANM to the equations ([l7| ) and having in mind the 
decompositions (^5[) , (|26|) , we attain the vector equation 



Au p + Av + Aw (t> = - (Ay + As) . 



(31) 



which represents an algebraic system consisting of three linear scalar equations with respect to the 
three unknowns p, lo and <j). 

The general sequence of the algorithm can be recapitulated in the following way. Let us assume 
that the functions yi.fe(a:), ye,fc(a;), and parameters Rs,k, 'Ps,k are given for fc > 0. We solve the 
linear BVPs (p7|), ( psj ) and thus we compute the functions 8^.^(2;), Ui^k{x), Vj_fe(a:), Wi^k{x) in the inner 
domain x g [0, 1]. Then we solve the linear BVPs (p9|), ( |30| ) in the external domain x G [l,oo] and 
compute the functions Se,kix),Ue^h{x), Ve,kix) and We.k{x)- Next, to obtain the increments pk, (^k 
and we solve the linear algebraic system (|3l]). So using the decompositions (p5|), ( p6| ) and then the 
formulae ( p^ - ( p2| ) we calculate the functions yi^fc+i(a;), ye,fc+i(x), the radius of the star Rs^k+i, the 
quantity ^k+i and the parameter boundary condition ips.k+i as well at the new iteration stage fc + 1. 

At the every iteration fc an optimal time step Topf is determinated in accordance with the Er- 
makov&Kalitkin formula [nO| 



S{0) 



<5(0) + <5(l) 

where the residual S{t) is calculated as follows 

S{Tk) = max [Sf, {Rs,k + TkPkT , i^k + TkUJkf , (VsM + Tk(f>kf 



and 



and Sf is the Euclidean residual of right hand side of the first equations in the systems (27), 

(§), (igl). 

The criterion for termination of the iterations is 6{Topt) < £ where e '--^ 10 ^ 10 for some fc. 

Taking into account the smoothness of sought solutions we solve the linear BVPs (p7|), ( p8| ) and 
(p9| ) , (^0|) employing Hermitean splines and spline collocation scheme of fourth order of approximation 
pi| . At that we utilize essentially the important feature that everyone of the above mentioned two 
groups vector BVPs (inner and external) have one and the same left hand sides. 

It is worth to note that the algebraic systems of linear equations and the system (^TJ) as well become 
ill-posed in the vicinity of the "exact" solution, i.e., for sufficiently small residuals S. That is why for 
small 6, for example ii 5 < 10^'^ (then Topt ^ 1 usually), it is expedient to use the Newton-Kantorovich 
method when the respective matrices are fixed for some 6 > 10""^. 



4 Some numerical results 



For a purpose of illustrating we will consider and discuss some results obtained from numerical exper- 
iments. A detailed description and analysing of results from physical point of view will be object of 
another our paper. 

In present article we consider concrete scalar-tensor model with functions (see Section ||) 
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Figure 1: The function (j{x) for ctc — 0.4; /ic = 1.2. 
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Figure 2: The function (p{x) for CTc = 0.4; /Zc = 1.2. 

The quantities b, A are given parameters. For completeness we note that in the concrete case the 
functions /(/i) and ^(/i) represent in parametric form the equation of state of noninteracting neutron 
gas while the function M^(cr^) describes the boson field with quartic self-interaction. 

The calculated functions <t{x), ip{x), iy{x) and ij{x) are plotted correspondingly in Fig. |l|, |2[ ^ andjj 
for the values of the parameters 7 = 0.1, A = 10 and 6=1. The behaviour of the mentioned functions 
is typical for wider range of the parameters not only for those values presented in the figures. The 
function cr(x) decreases rapidly from its central value Cc = 0.4 (in the case under consideration) to zero, 
at that when dimensionless coordinate x > 6 the function does not exceed 10"'*. Similarly the function 
i'{x) has most large derivative for x G (0, 9) after that it approaches slowly zero at infinity like K For 
example when x « 9 the derivative v'{x) ~ 10~^, while for x > 27 v'{x) < 10^'^, i.e., the asymptotical 
behaviour of calculated grid function and its derivative agrees very well with the theoretical prediction 
(see 1^). The function (p{x) increases rapidly for x < 4, besides that it trends asymptotically to 
zero. Obviously the quantitative behaviour of (p{x) for central value CTc = 0.4 is determinated by the 

B _F 

dominance of the term T over the term T (see [^). At last the function //(x) is nontrivial in the inner 
domain x € [0, 1], i.e. inside the star. Here it varies monotonously and continuously from its central 
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Parameters: A = 10, Y = 0.1, h = I 
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Figure 3: The function for CTc = 0.4; fic — 1.2. 



value (in the case under consideration) /ic = 1-2 till zero at a; = 1 corresponding to the radius of the 
star. 




Relative Coordinate jc 



Figure 4: The function for ac ~ 0.4; = 1-2. 



From physical point of view it is important to know the mass of the boson-fermion star and the 
total number of particles (bosons and fermions) making up the star. 
The dimensionless star mass can be calculated via the formula 

M = f -" {to" + To'+ exp(-A) J) ' + (^)^ dr. 

The dimensionless rest mass of the bosons (total number of bosons times the boson mass) is given 

by 

Mrb ^nj^ r^A\^) exp (^y^) '^^dr. 
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The dimensionless rest mass of the fermions is correspondingly 



Mrf = & / r^A'^iip) exp 
/o 



n(fi)dr 



where n{^) is the density of the fermions. In the case we consider we have n{^) = n'^{x). 

The dependencies of the star mass M (sohd hne) and the rest mass of fermions Mjip (dash hne) 
on the central value of the function are shown in configuration diagram on Fig. I for A = 0, 
7 = 0.1, b = 1 and CTc = 0.002. It should be pointed that for so small central value CTc we have 
in practice pure fermionic star. On the figure it is seen that from small values of fic to values near 
beyond the peak the rest mass is greater than the total mass of the star which means that the star is 
potentially stable. 
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- Star mass 

- Rest mass of fermions 



1 2 3 4 5 
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Figure 5: The star mass M and the rest fermion mass Mfip as functions of the central value fj,c 
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Figure 6: The binding energy E versus the rest fermion mass Mjip 

On Fig. ^ the binding energy of the star Eh = M — Mrb ~ Mj^p is drawn against the rest mass of 
fermions Mrp for A = 0, 7 = 0.1, 6 = 1 and ctc = 0.002. Fig. ^ is actually a bifurcation diagram. With 
increasing the central value of function one meets a cusp. The appearance of a cusp shows that 
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the stability of the star changes - one perturbation mode develops instability. Beyond the cusp the 
star is unstable and may collapse eventually forming black hole. The corresponding physical results 
for pure boson stars are considered in our recent paper 
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